!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   FEMilaro -- Finite Element Method toolkit
!!
!! FEMilaro is free software; you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 3 of the License, or
!! (at your option) any later version.
!!
!! FEMilaro is distributed in the hope that it will be useful, but
!! WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
!! General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with FEMilaro; If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>

!>\brief
!! 
!! Turbulent flux, Smagorinsky-Ricahrdson model (scalar viscosity)
!!
!! \n
!!
!! In this module we consider the model discussed in <a
!! href="http://onlinelibrary.wiley.com/doi/10.1111/j.2153-3490.1962.tb00128.x/pdf">
!! [Lilly, 1962]</a> and <a
!! href="http://journals.ametsoc.org/doi/pdf/10.1175/1520-0493%281983%29111%3C2341%3AACMFTS%3E2.0.CO%3B2">
!! [Durran, Klemp, 1983]</a>. In practice, this amounts to a
!! Smagorinsky eddy viscosity corrected according to the local
!! Richardson number.
!!
!! We set
!! \f{displaymath}{
!!  \underline{\underline{\mathcal{F}}}^d_{\underline{U}} =
!!  -\rho\left(\nu+\nu_{sgs}\right)\mathcal{S}, \qquad
!!  \underline{\mathcal{F}}^d_E = -c_p\rho\nu_E \nabla T.
!! \f}
!! with
!! \f{displaymath}{
!!  \mathcal{S} = \nabla\underline{u} + \nabla\underline{u}^T -
!!  \frac{2}{d}\nabla\cdot\underline{u}\,\mathcal{I}.
!! \f}
!! The momentum turbulent viscosity is
!! \f{displaymath}{
!!  \nu_{sgs} = f_{Ri}f_{y} \, C_s^2 \Delta^2 \,
!!  \frac{|\mathcal{S}|}{\sqrt{2}},
!! \f}
!! where \f$\Delta=\sqrt[3]{|K|}\f$, \f$
!! |\mathcal{S}| = \sqrt{\sum_{i,j}\mathcal{S}^2_{ij}}\f$,
!! and the two correction factors \f$f_{Ri}\f$ and \f$f_{y}\f$ account
!! for convective instabilities and wall effects, respectively.
!! Defining the Richardson number
!! \f{displaymath}{
!!  Ri = 2 \frac{\frac{g}{\theta} \frac{\partial\theta}{\partial
!!  z}}{|S|^2},
!! \f}
!! the first correction factor is
!! \f{displaymath}{
!!  f_{Ri} = \sqrt{\max\left(1-\frac{Ri}{Pr}\,,\sigma_{min}\right)},
!! \f}
!! where the coefficient \f$\sigma_{min}\f$ should be zero, which
!! would imply the possibility of completely suppressing turbulent
!! mixing; a value larger than zero is often useful for numerical
!! stability.

!! \f{displaymath}{
!!  \nu_U = k^2 \Delta x\,\Delta z\, \frac{|\mathcal{S}|}{\sqrt{2}}
!!  \sqrt{\max\left(1-\frac{Ri}{Pr}\,,\sigma_{min}\right)}, \qquad
!!  \nu_E = Pr^{-1}\nu_U,
!! \f}
!! where \f$k=0.21\f$, \f$Pr=\frac{1}{3}\f$ and
!! \f{displaymath}{
!!  |\mathcal{S}| = \sqrt{\sum_{i,j}\mathcal{S}^2_{ij}} , \qquad Ri =
!!  2 \frac{\frac{g}{\theta} \frac{\partial\theta}{\partial
!!  z}}{|S|^2}.
!! \f}
!! The coefficient \f$\sigma_{min}\f$ should be zero, which would
!! imply the possibility of completely suppressing turbulent mixing.
!! A value larger than zero is often useful for numerical stability.

!! The dissipative fluxes are obtained adding a viscous contribution,
!! computed as in \c mod_viscous_flux, and a turbulent contribution.
!! For the turbulent contribution, we use a dynamic model with
!! coefficients averaged on each element.
!!
!! Concerning the momentum flux, we have
!! \f{displaymath}{
!!  \underline{\underline{\mathcal{F}}}^d_{\underline{U}} =
!!  \underbrace{
!!  \underbrace{-\rho\nu_{sgs}\left[ \nabla\underline{u} + \nabla\underline{u}^T -
!!  \frac{2}{d}\nabla\cdot\underline{u}\,\mathcal{I}\right]}_{
!!  \displaystyle\underline{ \underline{\tau}}^d}
!!  +\underbrace{C_I\rho\Delta^2 \left|\nabla\underline{u} +
!!  \nabla\underline{u}^T\right|^2\,\mathcal{I}}_{
!!  \displaystyle\frac{1}{d}\rho e_K\mathcal{I}} }_{\displaystyle 
!!  \displaystyle\underline{ \underline{\tau}}}
!!  +\,\underline{\underline{\mathcal{F}}}^{visc}_{\underline{U}} 
!! \f}
!! where \f$d\f$ is the spatial dimension, \f$\Delta\f$ is a
!! characteristic scale,
!! \f{displaymath}{
!! \nu_{sgs} = C_D\Delta^2 \left|\nabla\underline{u} +
!!  \nabla\underline{u}^T\right|,
!! \f}
!! and the viscous terms are described in \c mod_viscous_flux. Notice
!! that this choice ensures \f${\rm
!! tr}(\underline{\underline{\tau}}^d)=0\f$, so that the trace of the
!! turbulent stress tensor is \f$\rho e_K\f$, which can be regarded as
!! the subgrid kinetic energy. The two dimensionless constants
!! \f$C_D\f$ and \f$C_I\f$ are determined dynamically. To this
!! end, let us introduce the following notation:
!!
!<----------------------------------------------------------------------
module mod_smagrich
!General comments:
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------

 use mod_used, only: &

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_smagrich_constructor, &
   mod_smagrich_destructor,  &
   mod_smagrich_initialized

 private

!-----------------------------------------------------------------------

! Module types and parameters

 ! public members
 logical, parameter ::             &
 integer, parameter ::             &
 ! private members
 real, parameter ::                &

! Module variables

 ! public members
 logical, protected ::               &
   mod_smagrich_initialized = .false.
 integer ::               &
 ! private members
 real ::                  &

 character(len=*), parameter :: &
   this_mod_name = 'mod_smagrich'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_smagrich_constructor()
  integer, intent(in) :: &
  real, intent(in) ::    &

  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_used1.eqv..false.) .or. &
       (mod_used2.eqv..false.) ) then
   endif
   if(mod_smagrich_initialized.eqv..true.) then
   endif
   !----------------------------------------------

   mod_smagrich_initialized = .true.
 end subroutine mod_smagrich_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_smagrich_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_smagrich_initialized.eqv..false.) then
   endif
   !----------------------------------------------

   mod_smagrich_initialized = .false.
 end subroutine mod_smagrich_destructor

!-----------------------------------------------------------------------

end module mod_smagrich

